#ifndef UHM_SOLVER_BASE_SOLVER_HXX_
#define UHM_SOLVER_BASE_SOLVER_HXX_

#define UHM_TREE_MAX_CHILDREN 3

namespace uhm {
  // Forward Decalaration
  template<class S_> 
  Int_ op_leaf_to_me_seq(S_ *s, typename S_::_Tree_::Tnd_ *me,
                         Int_ (*op_func)(S_*, typename S_::_Tree_::Tnd_*));
  template<class S_> 
  Int_ op_me_to_leaf_seq(S_ *s, typename S_::_Tree_::Tnd_ *me,
                         Int_ (*op_func)(S_*, typename S_::_Tree_::Tnd_*));
  template<class S_> 
  Int_ op_leaf_to_me_par(S_ *s, typename S_::_Tree_::Tnd_ *me,
                         Int_ (*op_func)(S_*, typename S_::_Tree_::Tnd_*));
  template<class S_> 
  Int_ op_me_to_leaf_par(S_ *s, typename S_::_Tree_::Tnd_ *me,
                         Int_ (*op_func)(S_*, typename S_::_Tree_::Tnd_*));

  class Solver_Base_ : public Disp_ {
  public:

  protected:
    Int_ _type, _uplo, _sym, _rhs, _blk, _caching, _offset, _level;
    
  public:
    Solver_Base_() : _type(Mat_Base_::DOUBLE),
                     _uplo(Mat_Base_::FULL),
                     _sym(Mat_Base_::UNSYM),
                     _rhs(1),
                     _blk(256),
                     _caching(0),
                     _offset(0),
                     _level(-1) { }

    Solver_Base_(Int_ type, Int_ uplo, Int_ sym,
                 Int_ rhs, Int_ blk) :  _type(type),
                                        _uplo(uplo),
                                        _sym(sym),
                                        _rhs(rhs),
                                        _blk(blk),
                                        _caching(0),
                                        _offset(0),
                                        _level(-1) { }
    
    Solver_Base_(Int_ type, Int_ rhs, Int_ blk) :  _type(type),
                                                   _uplo(Mat_Base_::FULL),
                                                   _sym(Mat_Base_::UNSYM),
                                                   _rhs(rhs),
                                                   _blk(blk),
                                                   _caching(0),
                                                   _offset(0),
                                                   _level(-1) { }
    virtual~Solver_Base_() { }

    virtual std::ostream& content(std::ostream& os) const  {
      return (os << "- Solver Basic Info -\n" 
              << "Type(" << _type << "), "
              << "Uplo(" << _uplo << "), "
              << "Sym(" << _sym << "), "
              << "Rhs(" << _rhs << "), "
              << "Blk(" << _blk << ")\n");
    }
    virtual std::ostream& report(std::ostream& os)   {
      return (os << "- Report - \n" << "Nothing to Report\n");
    }

    void set_offset(Int_ offset) { _offset = offset; }
    Int_ get_offset() const { return _offset; }

    void set_level(Int_ level) { _level = level; }
    Int_ get_level() const { return _level; }

    Int_ get_type() const { return _type; }
    Int_ get_uplo() const { return _uplo; }
    Int_ get_sym() const { return _sym; }
    
    Int_ get_rhs() const { return _rhs; }
    Int_ get_blk() const { return _blk; }

    void set_caching(Int_ caching) { _caching = caching; }
    Int_ get_caching() const { return _caching; }

    // Pure ...
    virtual Int_ copy_in_begin()  { return false; }
    virtual Int_ copy_in_end()  { return false; }

    virtual Int_ copy_in(Int_ *el, Int_ eids,
                         Int_ *nodes, Int_ nids, 
                         Int_ *dofs, Int_ *weights, Int_ *kinds, Int_ n_nodes,
                         Int_ cs, void *buf_A, void *buf_B)  { return false; }
    virtual Int_ copy_in(Int_ *nodes, Int_ ids, 
                         Int_ *dofs, Int_ n_nodes, 
                         Int_ cs, void *buf_B)  { return false; }
    virtual Int_ copy_out(Int_ *el, Int_ ids, Int_ cs, void *buf_X)  { return false; }
    virtual Int_ copy_out(Int_ *nodes, Int_ ids, Int_ n_nodes, 
                          Int_ cs, void *buf_X)  { return false; }
    virtual Int_ copy_out_root(Int_ *nodes, Int_ ids, Int_ n_nodes, 
                               Int_ cs, void *ABR)  { return false; }
    virtual Int_ create()  { return false; }
    virtual Int_ clear()  { return false; }
    virtual Int_ flush()  {  return false; }
  };

  template<class Id_Nd_, class Id_El_, class Uhm_>
  class Solver_ : public Solver_Base_ {
  public:    
    typedef Solver_<Id_Nd_,Id_El_,Uhm_> _Solver_;
    typedef Mesh_<Id_Nd_,Id_El_,Uhm_> _Mesh_;
    typedef Tree_<typename _Mesh_::Element_*,UHM_TREE_MAX_CHILDREN> _Tree_;
    typedef Cache_<Id_Nd_,Uhm_> _Cache_;

    class Stat_ : public Disp_ {
    private:
      Double_ _flop, _recorded, _storage;
      Int_ _dof, _min, _max, _n;
    public:
      Stat_() : _flop(0.0), _recorded(0.0), _storage(0.0), 
                _dof(0), _min(0), _max(0), _n(0) { }
      virtual~Stat_() { }
      void update(Stat_ &b) {
        b._flop += _flop; b._recorded += _recorded; b._storage += _storage;
        b._dof += _dof;
        b._n += _n;

        if (b._min)
          b._min = (b._min > _min ? _min : b._min);
        else
          b._min = _min;
        if (b._max)
          b._max = (b._max < _max ? _max : b._max);
        else
          b._max = _max;
      }
      void update(typename _Mesh_::Element_ *e) {
        Uhm_ &uhm = e->uhm(); 
        
        _flop     += uhm.get_leading_flop(); 
        _recorded += e->get_flop();

        _storage  += uhm.get_allocated_storage();
        _dof      += uhm.get_mt(); 
        
        Int_ mm = (uhm.get_mt() + uhm.get_mb());

        if (_min == 0) _min = mm;
        _min = (_min > mm ? mm : _min); 

        if (_max == 0) _max = mm;
        _max = (_max < mm ? mm : _max);

        ++_n;
      }
      void clear() {
        _flop = 0.0; _recorded = 0.0; _storage = 0.0;
        _dof  = 0; _min = 0; _max = 0;
        ++_n;
      }
      virtual std::ostream& content(std::ostream& os) const  {
        return (os 
                << "  FLOP    = " << (_flop/1.0e9) << ", "
                << "  Updated = " << (_recorded/1.0e9) << " [GFLOP]\n"
                << "  Storage = " << (_storage/1.0e6)  << " [MB]\n"
                << "  DOF     = " << _dof << ","
                << "  # UHMs  = " << _n << ","
                << "  Min UHM = " << _min << ","
                << "  Max UHM = " << _max << "\n");
      }
    };

  protected:
    _Mesh_  _mesh;
    _Tree_  _tree;
    _Cache_ _work;
    _Cache_ _cache;

    std::vector<typename _Tree_::Tnd_*> _tnds;
    Map_<Id_Nd_> _nodes;

    std::vector<Stat_> _stat;
    
  public:
    Solver_() : Solver_Base_() {  } 
    Solver_(Int_ type, Int_ uplo, Int_ sym,
            Int_ rhs, Int_ blk) :  Solver_Base_(type, uplo, sym, rhs, blk) { }
    Solver_(Int_ type, Int_ rhs, Int_ blk) : Solver_Base_(type, rhs, blk) { } 
    virtual~Solver_() { }

    virtual std::ostream& content(std::ostream& os) const  {
      this->Solver_Base_::content(os);
      os << "\n";
      os << "  Work\n"
         << "    # of items    = " << _work.get_size() << "\n"
         << "  Cache\n"
         << "    # of items    = " << _cache.get_size() << "\n" 
         << "    caching level = " << this->get_caching() << "\n";

      os << "  Tree\n"
         << "    # of items    = " << _tree.get_n_nodes() << "\n";

      if (_tree.get_root() != NULL)
        os << "    height        = " << _tree.get_root()->get_gen() << "\n";

      os << "  Mesh\n"
         << "    # of nodes    = " << _mesh.get_n_nodes() << "\n"
         << "    # of elements = " << _mesh.get_n_elements() << "\n";

      if (_tnds.size()) 
        os << "  Active TreeNodes\n"
           << "    # of items    = " << _tnds.size() 
           << "    cutting level = " << this->get_level() << "\n";
      return os;
    }

    virtual std::ostream& report(std::ostream& os)   {
      if (_tree.get_root() == NULL || !_tree.is_sorted()) 
        return (os << "Nothing to report\n");
      
      // Show basic info of this solver
      this->content(os);
      _stat.clear();

      Int_ size = (this->get_level() < 0 ? _tree.get_root()->get_gen() : this->get_level());
      ++size;
      _stat.reserve(size);
      for (Int_ i=0;i<size;++i)
        _stat.push_back(Stat_());
      
      for (Int_ i=0;i<_tree.get_size();++i) {
        typename _Tree_::Tnd_ *tnd = _tree.get_tnd_from_sorted(i);
        if (tnd->get_gen() < size)
          _stat[tnd->get_gen()].update(tnd->get_data());
      }

      Stat_ tot;
      for (Int_ i=0;i<size;++i) {
        os << "\n\n" << i << " Level\n" << _stat[i];
        _stat[i].update(tot);
      }

      os << "------------------------------------\n";
      os << "  Total\n" << tot << "\n";

      Double_ storage, max_storage;
      storage_counter(storage, max_storage);

      os << "  Storage\n"
         << "    Work     = " << (_work.get_allocated_storage()/1.0e6) << " [MB] \n"
         << "    Cache    = " << (_cache.get_allocated_storage()/1.0e6) << " [MB] \n\n" 
         << "    Current  = " << (storage/1.0e6) << " [MB] \n" 
         << "    Max Used = " << (max_storage/1.0e6) << " [MB] \n";

      return (os << "\n");
    }

    virtual void get_flop(Double_ &flop, Double_ &recorded) {
      UHM_ERROR((_tree.get_root() != NULL && _tree.is_sorted()), UHM_ERROR_SOLVER_GET_FLOP);
      flop = 0.0; recorded = 0.0;
      for (Int_ i=0;i<_tree.get_size();++i) {
        typename _Tree_::Tnd_ *tnd = _tree.get_tnd_from_sorted(i);
        if (this->get_level() < 0 || tnd->get_gen() <= this->get_level()) {
          flop     += tnd->get_data()->uhm().get_leading_flop();
          recorded += tnd->get_data()->get_flop();
        }
      }
    }

    // ** Basic interfaces
    _Mesh_&  mesh()  { return _mesh; }
    _Tree_&  tree()  { return _tree; }
    _Cache_& work()  { return _work; }
    _Cache_& cache() { return _cache; }

    // I do not want to open access to internal variables, but
    // for experimental puspose... now they are accessable.
    std::vector<typename _Tree_::Tnd_*>& tnds() { return _tnds; }
    Map_<Id_Nd_>& nodes() { return _nodes; }
    
    // ** Application mapping interfaces
    virtual void copy_get_ready(std::vector< std::pair<Id_Nd_,Int_> > &nodes,
                                Map_<Id_Nd_> &factor, Map_<Id_Nd_> &schur) {
      factor.push_begin(0);
      factor.push_end();

      schur.push_begin(nodes.size());
      for (Int_ i=0;i<nodes.size();++i) {
        schur.push(nodes[i].first);
      }
      schur.push_end(); 
    }

    // ** C++ Interfaces
    using Solver_Base_::copy_in;

    virtual Int_ copy_in(Id_El_ el, 
                         std::vector< std::pair<Id_Nd_,Int_> > &nodes, 
                         Int_ cs, void *buf_A, void *buf_B) {
      // Add element to mesh, nodes are stored and sorted 
      typename _Mesh_::Element_* e = _mesh.add_element(el);
      e->push_begin();
      for (Int_ i=0;i<nodes.size();++i) {
        typename _Mesh_::Node_* n = _mesh.add_node(nodes[i].first, nodes[i].second);
        e->push(n);
      }
      e->push_end();

      // Copy matrices to work, do not use nds2map
      // nds in element is already sorted.
      Map_<Id_Nd_> factor, schur; copy_get_ready(nodes, factor, schur); 
      Int_ dof = schur.dof();

      // dummy is always required when we use push, pull is okay without dummy
      Uhm_ dummy(this->_type, this->_uplo, this->_sym, 0, dof, this->_rhs, this->_blk);
      typename _Cache_::Entry_ in(factor, schur, dummy), *out;
      UHM_ERROR(_work.push(in, &out), UHM_ERROR_CACHE_PUSH);

      // Create matrices at a schur storage
      out->uhm().create_flat_without_buffer();
      out->uhm().create_flat_buffer(true);

      typename Uhm_::_A_ A(this->_type, dof, dof);  
      typename Uhm_::_B_ B(this->_type, dof, _rhs);

      A.create_flat_without_buffer();         B.create_flat_without_buffer();
      A.attach_flat_buffer(cs, 1, buf_A);     B.attach_flat_buffer(cs, 1, buf_B);
      copy(A, out->uhm().ABR());              copy(B, out->uhm().BB());
      A.free_flat_without_buffer();           B.free_flat_without_buffer();

      return true;
    }

    // Update RHS only
    virtual Int_ copy_in(std::vector< std::pair<Id_Nd_,Int_> > &nodes, 
                         Int_ cs, void *buf_B) {
      // Copy matrices to work
      Map_<Id_Nd_> factor, schur; copy_get_ready(nodes, factor, schur);
      Int_ dof = schur.dof();

      typename _Cache_::Entry_ in(factor, schur), *out; 
      UHM_ERROR(_work.pull(in, &out), UHM_ERROR_SOLVER_PULLED_WORK); 

      typename Uhm_::_B_ B(this->_type, dof, this->_rhs);

      B.create_flat_without_buffer();
      B.attach_flat_buffer(cs, 1, buf_B);
      copy(B, out->uhm().BB());
      B.free_flat_without_buffer();

      return true;
    }
    
    virtual Int_ copy_out(std::vector< std::pair<Id_Nd_,Int_> > &nodes,
                          Int_ cs, void *buf_X) {
      Map_<Id_Nd_> factor, schur; copy_get_ready(nodes, factor, schur);
      Int_ dof = schur.dof();

      typename _Cache_::Entry_ in(factor, schur), *out;
      UHM_ERROR(_work.pull(in, &out), UHM_ERROR_SOLVER_PULLED_WORK);

      typename Uhm_::_B_ X(this->_type, dof, this->_rhs);
      X.create_flat_without_buffer();
      X.attach_flat_buffer(cs, 1, buf_X);
      copy(out->uhm().XB(), X);
      X.free_flat_without_buffer();

      return true;
    }
    
    virtual Int_ copy_out(Id_El_ el, Int_ cs, void *buf_X) {
      typename _Mesh_::Element_ *e = _mesh.find_element(el);
      Map_<Id_Nd_> factor, schur;
      e->nds2map(factor, schur);
      Int_ dof = schur.dof();

      typename _Cache_::Entry_ in(factor, schur), *out;
      UHM_ERROR(_work.pull(in, &out), UHM_ERROR_SOLVER_PULLED_WORK);

      typename Uhm_::_B_ X(this->_type, dof, this->_rhs);
      X.create_flat_without_buffer();
      X.attach_flat_buffer(cs, 1, buf_X);
      copy(out->uhm().XB(), X);
      X.free_flat_without_buffer();

      return true;
    }
    
    virtual Int_ copy_out_root(std::vector< std::pair<Id_Nd_,Int_> > &nodes,
                               Int_ cs, void *ABR) {
      return false;
    }


    // Generic interfaces
    virtual Int_ flush()     { 
      _work.clear(); 
      _cache.clear();
      
      return true; 
    }

    virtual Int_ copy_in_begin() { 
      _tree.clear();
      _work.clear();
      _mesh.clear();

      time_in("uhm:: solver:: copy_in (begin - end)", UHM_TIME_LEVEL_FRONT);

      _work.push_begin();

      return true;
    }

    virtual Int_ copy_in_end() {
      _work.push_end();

      time_out(UHM_TIME_LEVEL_FRONT);

      return true;
    }

    virtual Int_ copy_in(Int_ *elt, Int_ eids,
                         Int_ *nodes, Int_ nids, 
                         Int_ *dofs, Int_ *weights, Int_ *kinds, Int_ n_nodes,
                         Int_ cs, void *buf_A, void *buf_B) {
      Id_El_ el;
      for (Int_ i=0;i<eids;++i) 
        el.set_id(i, elt[i]);

      std::vector< std::pair<Id_Nd_,Int_> > nodes_local; 
      nodes_local.reserve(n_nodes);

      Id_Nd_ nd;
      for (Int_ j=0;j<n_nodes;++j) {
        for (Int_ i=0;i<nids;++i) 
          nd.set_id(i, nodes[j*nids+i]);
        nd.set_dof(dofs[j]);
        nd.set_weight(weights[j]);
        nodes_local.push_back(std::make_pair(nd, kinds[j]));
      }

      return copy_in(el, nodes_local, cs, buf_A, buf_B);
    }

    virtual Int_ copy_in(Int_ *nodes, Int_ ids,
                         Int_ *dofs, Int_ n_nodes,
                         Int_ cs, void *buf_B) {
      std::vector< std::pair<Id_Nd_,Int_> > nodes_local;
      nodes_local.reserve(n_nodes);

      Id_Nd_ nd;
      for (Int_ j=0;j<n_nodes;++j) {
        for (Int_ i=0;i<ids;++i)
          nd.set_id(i, nodes[j*ids+i]);
        nd.set_dof(dofs[j]);
        nodes_local.push_back(std::make_pair(nd, 0));
      }

      return copy_in(nodes_local, cs, buf_B);
    }
    virtual Int_ copy_out(Int_ *elt, Int_ ids, Int_ cs, void *buf_X) {
      Id_El_ el;
      for (Int_ i=0;i<ids;++i)
        el.set_id(i, elt[i]);

      return copy_out(el, cs, buf_X);
    }
    virtual Int_ copy_out(Int_ *nodes, Int_ ids, Int_ n_nodes,
                          Int_ cs, void *buf_X) {
      std::vector< std::pair<Id_Nd_,Int_> > nodes_local;
      nodes_local.reserve(n_nodes);

      Id_Nd_ nd;
      for (Int_ j=0;j<n_nodes;++j) {
        for (Int_ i=0;i<ids;++i)
          nd.set_id(i, nodes[j*ids+i]);
        nodes_local.push_back(std::make_pair(nd, 0));
      }

      return copy_out(nodes_local, cs, buf_X);
    }
    virtual Int_ copy_out_root(Int_ *nodes, Int_ ids, Int_ n_nodes,
                               Int_ cs, void *ABR) { 
      return false;
    }

    // Create tree for mesh. 
    // Matrices in work are already fullfilled by copying operations.
    virtual Int_ create()    { 
      time_in("uhm:: solver:: create", UHM_TIME_LEVEL_FRONT);
      time_in("uhm:: solver:: create:: Mesh2Tree",UHM_TIME_LEVEL_INTERNAL);
      mesh2tree<_Mesh_,_Tree_>(this->mesh(),this->tree());
      time_out(UHM_TIME_LEVEL_INTERNAL);

      time_in("uhm:: solver:: create:: Separating Nodes and Creating Matrices",UHM_TIME_LEVEL_INTERNAL);
#pragma omp parallel
      {
#pragma omp single nowait
        {
          UHM_ERROR(op_tree_level_par<_Tree_>(_tree, &op_elt_separate_nodes<_Tree_>),
                    UHM_ERROR_SOLVER_SEPARATE_NODES);
          UHM_ERROR(op_me_to_leaf_seq<_Solver_>(this, _tree.get_root(),
                                                &op_uhm_set_offset<_Solver_>),
                    UHM_ERROR_SOLVER_SET_OFFSET);

          if (_level >= 0) {
            _tree.get_tnds_from_sorted(_level, _tnds);
          } else {
            _tnds.push_back(_tree.get_root());
          }

          for (Int_ i=0;i<_tnds.size();++i) {
#pragma omp task firstprivate(i)
            {
              UHM_ERROR(op_leaf_to_me_par<_Solver_>(this, _tnds[i],
                                                    &op_uhm_set_and_create_matrices<_Solver_>),
                        UHM_ERROR_SOLVER_CREATE_MATRICES);
            }
          }
        }
      }
      time_out(UHM_TIME_LEVEL_INTERNAL);

      time_out(UHM_TIME_LEVEL_FRONT);
      return true;
    }

    // Clear all object except for cache. 
    // In adaptive sequence, this enables to prepare for a next iteration.
    virtual Int_ clear()     { 

      time_in("uhm:: solver:: clear",UHM_TIME_LEVEL_FRONT);
      _tree.clear();
      _mesh.clear();
      _tnds.clear();
      time_out(UHM_TIME_LEVEL_FRONT);

      return true; 
    }
  };
}


#endif
